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SUMMARY 

We consider the applications of stress-strength models in studies involving clinical 
trials. When studying the effects and side effects of certain procedures (treatments), it 
is often the case that observations are correlated due to subject effect, repeated measure- 
ments and observing many characteristics simultaneously. We develop maximum likelihood 
estimator (MLE) and uniform minimum variance unbiased estimator (UMVUE) of the re- 
liability which in clinical trial studies could be considered as the chances of increased side 
effects due to a particular procedure compared to another. The results developed apply to 
both univariate and multivariate situations. Also, for the univariate situations we develop 
simple to use lower confidence bounds for the reliability. Further, we consider the cases 
when both stress and strength constitute time dependent processes. We define the future 
reliability and obtain methods of constructing lower confidence bounds for this reliability. 
Finally, we conduct simulation studies to evaluate all the procedures developed and also 
to compare the MLE and the UMVUE. 
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1. Introduction: 

Critically ill patients often require endotracheal intubation and mechanical ventilation. 
Endotracheal intubation is a procedure by which a tube that is attached to a ventilation 
machine is inserted into the trachea through the mouth of the patient, allowing oxygen to 
pass through the airway. The presence of the tube interferes with the mucociliary transport 
system - the normal mechanism by which respiratory tract secretions are removed (see 
Gray et al. (1990)). Thus, endotracheal suctioning is an important component of the 
mechanically ventilated patient’s care. In the presence of pulmonary disease, infection 
or dehydration, respiratory secretions may become thick and tenacious, making removal 
by suctioning difficult. Therefore, instillation of normal saline solution (5-10 ml bolus 
of sterile-normal-saline) into the tube prior to suctioning is a common practice. This 
procedure is believed to aid in eliciting a cough and also liquefy and mobilize secretions. 
Several studies have been conducted to study the side effects of normal saline instillation 
(NSI) on variables such as gas exchange, respiratory system, etc. However, Gray et al. 
(1990) were the first to conduct a designed clinical study. Further, Jordan and Garrett 
(1991) conducted another clinical trial to study the effect of NSI on variables such as heart 
rate, ventilatory rate and Sa0 2 (the amount of oxygen content in hemoglobin). The effect 
of NSI can be studied in two ways: (i) compare these means of the variables with NSI and 
without NSI or (ii) estimate the chances that a patient gets worse (or better), in terms of 
these variables, by using NSI compared to not using NSI. Here, we are interested in the 
latter. For example, we may want to assess the chances of increased heart rate due to NSI 
procedure compared to without NSI procedure. In other words, if yu denotes the increase 
in the heart rate of the ith subject with NSI procedures, and y 2i denotes that without the 
NSI procedures, then we are interested in 


R = P[Y u >Y 2i }. (1) 

This probability is usually known as the “reliability” in reliability stress-strength 
literature (see Johnson (1988), Guttman et al. (1988), Weerahandi and Johnson (1992) 
and the references therein). One may compare the means but as noted by Weerahandi 
and Johnson (1992) ‘even when comparing two treatments, it may be more informative to 
study the unit free quantity, R, rather than comparing the means’. Especially when yu 
and y 2 i are clinical variables that help measure risk of a particular procedure, R would be 
of much more interest. 

Much of the work in the literature has been concentrated around situations where 
stress and strength are independent, univariate normal random variables (also see Church 
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and Harris (1970), Downtown (1971), Enis and Geisser (1971)). Owen, Craswell and Han- 
son (1964) developed upper nonparametric confidence bounds for the reliability by relaxing 
the distributional assumptions. Weerahandi and Johnson (1992) provided procedures for 
conducting a test such as H 0 : R < R 0 . In the presence of explanatory variables, methods 
of constructing confidence intervals for R have been discussed in Guttman et al. (1988). 
Most of the applications considered in the literature have been from the industrial point 
of view. We feel that these reliability stress-strength results can also be applied to the 
clinical data just as easily. However, in clinical setting, it is often the case that stress and 
strength would be correlated because of the subject effect, repeated measurements and 
observing many characteristics simultaneously. For the cases when stress and strength 
follow a bivariate normal distribution, Owen et al. (1964) provided confidence intervals for 
the reliability via the non-central t-distribution. This of course includes the paired data 
sets, where stress and strength are the paired observations on various subjects. 

In this paper we first consider the case of bivariate normal distribution (as in Owen 
et al. (1964)), and provide the maximum likelihood estimator (MLE) and also give a 
procedure of constructing UMVUE of R. Since MLE and UMVUE are not analytically 
comparable, we conduct simulation studies to compare their performance in terms of bias 
and mean square error (mse). It will be shown that UMVUE performs better when the true 
probability (i.e. reliability) is roughly between 0.25 and 0.75. We also give three easy to use 
approximations to construct (1 - a) 100% most accurate lower bounds for R. Simulation 
studies are used to evaluate and compare the approximations in terms of their coverage 
probabilities. Based on the evaluations, we recommend a particular approximation which 
maintains the coverage almost identical to the specified coverage. 

Next, we study the case when y and x are p-variate normal vectors observed on 
same patients using two different treatments or at two different time points using a single 
treatment (for example: before and after the treatment). Therefore, we allow y and x 
to be correlated. Once again, we provide the MLE and devise a procedure to obtain the 
UMVUE of R. Because of the multivariate nature of the problem, we do not get a closed 
form for the UMVUE. However, as long as the dimension of the vector is not too high, a 
simple program utilizing IMSL or NAG subroutines can be used to obtain the UMVUE. 
We compare MLE and UMVUE via simulation studies in terms of the bias and the mse. 
Simulation studies once again indicate that the UMVUE performs better than the MLE 
when the true probability is roughly between 0.25 and 0.75. Even though MLE has a closed 
form, and is easier to compute than the UMVUE, the advantages of UMVUE outweigh 
the complexity of its usage. For a good discussion on the advantages UMVUE in similar 
settings see Govindarajulu (1995). 
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Further, we deal with the situations when both stress and strength are observed 
simultaneously, at time points l,2,...,n, and they are time dependent. Assuming that 
the conditional joint distribution of (y n + i,®n+i), given all the past observations, depends 
only on the current observations y n and x n , we develop procedures to predict the reliability 
at time (n + 1). We also provide confidence intervals for the reliability at time (n + 1) 
conditional on the past observations. 

Finally, we apply the results of the first two procedures to the clinical trials data 
provided by Jordan and Garrett (1991). We also provide examples for the third procedure. 

2. Point and interval estimation of reliability, when stress and strength follow 
a bivariate normal distribution. 

Assume that y and x jointly follow a bivariate normal distribution given by, 

Q~JV(/i,£) (2) 

where p is the mean vector and E is the variance covariance matrix. 

We are interested in the reliability, defined by 

Ri = P[y > x] = P[y - x > 0]. (3) 

Let c T = (1,-1), then R = <f>((5), where $ denotes the standard normal distribution 
function and 5 = c'p/(c'Ec) 1/2 . Let ft and E denote the MLE’s of fi and E, based on the 
observations (yi, • • • ,y n ) and (xi, • • • ,x n ). Then, the maximum likelihood estimator of R 
is given by 

R f = <L(<5) (4) 

where 6 is <5 with y, and E replaced by their MLE’s, i.e 6 — P ///(c'Ec) 1 / 2 . 

Now we derive the UMVUE of R x . Let A denote an orthogonal matrix with its 
first two rows given by of = (1 — , ”) and aj = (^, • • • , ^)- Also, define 

I {y i_ X i) = 1 if ( yi - xi) > 0 and I( yi - Xl ) = 0 otherwise. Then, it is easy to see that 
hm~ x i) is an unbiased estimator of R\. Now, by making an orthogonal transformation 
y = Ax, and using Rao-Blackwell Theorem (for properties of complete sufficient statistics) 
and Basu’s Theorem (on ancillary statistics) (see Lehman (1983)) we obtain the UMVUE 
of R x given by 


y/n — 1 <5 

b umvue = I 1{z)dZ) 

— oo 


( 5 ) 
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where 


/(*) = 



<y¥ 

otherwise, 


with 6 1 = r(V)/( r (^) r (^))- 
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Lower confidence bounds for R\\ 

Using the standard distribution theory (see Rohatgi (1976)), it can be shown that 
c T /t ~ N(c T pi,±c T Zc), n (c T Ec) _1 (c r Ec) ~ xl-\ and hence, y/n - 1*5 ~ t n -i{y/nS), 
where t n -\{^Jn8) denotes the non-central ^-distribution with (n-1) d.f. and non-centrality 
parameter y/n8. If we use the unbiased estimator of E instead of the MLE, then y/n8 ~ 
tn- i(Vn6), where 6 is 8 with pi and E replaced by their unbiased estimators. 

Further, since constructing a lower bound on R\ is equivalent to constructing a lower 
bound on <5, we obtain a lower bound on 8. A uniformly most accurate, scale invariant, 
(1 — a) 100% lower confidence bound for 8 is given by the solution of 


F t n -i(JZ5){Vn8) - 1 - a, ( 6 ) 

where Ft . (■) denotes the distribution function of non-central ^-distribution with 

l n — 1 (v/SS) W 

(n — 1) d.f. and non-centrality parameter y/n8. Note that in (6) we have used the unbi- 
ased estimator of E instead of MLE. If one uses MLE, one needs to adjust the constant 
accordingly. 

Even though, solution of equation (6) is computationally tenable, a closed form ap- 
proximation that is simple to use and satisfies the necessary coverage properties is desirable 
in practice. Therefore, we provide three different closed form approximations and obtain 
the approximate lower bound in each case. 

Approximate lower confidence bounds: 

(i) Using a well known approximation of a non-central ^-distribution by a standard normal 
distribution (see Johnson and Kotz (1970)), we provide an approximate lower bound 
for 8 given by 


L Al (8) = 8b 2 



1 + 


n 

n — 1 



(7) 


where b 2 = b 3 = (n - 1)[1 - ^(T(f)/r( V)) 2 ]. “ d = 

$ -1 (l — a). If we further approximate this using \/2(r( ; |)/r( Ii ^d.)) = y'n — 1, our ap- 
proximation of non-central t by standard normal reduces to the one used in Guttman et 
al. (1988, Section 2). However, we will use the lower bound in (7) as it is not difficult to 
use. 


(ii) Another approximation is of the form 


— sinh -1 

0 4 



^5 


sinh 


-l 


'b 4 E{^8)\ 

h ) 


~N( 0,1), 
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and hence the lower bound is given by 


La,0) = 




^ 6 ) 


( 8 ) 


where 64 = r(^){2[(n - 3)r 2 (^) - l]" 1 }*, 65 = \A" “ x )/( n 3 ) and 

i>6 = 

(iii) A further approximation of the approximation in (ii) can be made to obtain a much 
simpler lower bound given by 


LaM = 


\fn 


\Jn — 1 sin h 


\/n8 
y/n — 1 


- Zr 


( 9 ) 


We should note that computation of the left hand side of the equation in (6) can be 
done via numerical integration without too much difficulty (as noted also by Guttman 
et al. (1988)). However, the approximations provide simple lower bounds that can be 
calculated easily and one only needs the standard normal distribution tables to compute 
the bounds. As it is clear that there is no direct analytical comparison between these 
three approximations, we conduct a simulation study and evaluate these approximations 
in terms of the coverage probability in the current setting. 


Simulation Studies: 

We generated random samples of sizes n=10, 15,20,25,30, and 35 from normal distribu- 
tion with various combinations of means and variances. For each simulation we computed 
the MLE, the UMVUE, and their respective biases. We also calculated the coverage prob- 
abilities using: equation in (6), approximation (i) , approximation (ii), and approximation 
(iii). At the end of 20,000 simulations, we find the averages and the standard deviations 
of these quantities. Further, we calculate the mean square error (mse) for both MLE and 
UMVUE. In Table I, we report the true probability, MLE, UMVUE. Biases, Standard 
deviations and the root mean square errors of MLE and UMVUE are also given in the 
same table. It is apparent from the results that UMVUE does better than MLE if the true 
probability is neither low or nor high (i.e. if true probability is roughly between 0.25 and 
0.75), otherwise MLE seems to do better. This is not very surprising because, for instance, 
when the true probability is very low (or high), even though MLE is biased, the bias will 
not contribute a great deal to the mse. Whereas, when the true probability is neither high 
nor low, the bias plays an important role and hence UMVUE performs better in terms 
of bias and mse. When an unbiased estimator is of interest, one should choose UMVUE 
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regardless of the true probability because when it does have more mse than that of the 
MLE, the mse is marginally more. However, if one prefers to use the estimate only based 
on smaller mse, it helps if one has some knowledge of the true probability (in terms of 
low, medium, high). In practice, one can calculate <5 and if this value is a high positive (or 
high negative) number, that indicates that the true probability is probably high (or low). 
This knowledge can be utilized in deciding between the two estimators. In clinical studies, 
usually, the cases that are of most interest are when the true probabilities are neither too 
high nor too low. As for the ease of use, certainly MLE is easier to calculate than the 
UMVUE, however the calculation of the XJMVXJE in this case is not at all difficult. 

INSERT TABLE I HERE 

Further, the true coverage (confidence level) and the simulated coverages obtained by 
the four methods mentioned above are reported in Table II. From this, it is also clear that 
the simulated coverages obtained using (6) and using approximation (i) are very close to 
each other and to the true coverage. Approximation (ii) seems to be very conservative 
in the sense of having lot more coverage than the desired confidence. In other words, the 
lower bound obtained using this approximation is much smaller than the true lower bound, 
resulting in a wider interval than necessary. Approximation (iii) does very badly until the 
true probability gets to be around 0.25. And as the true probability increases beyond 0.75, 
this approximation turns out to be very conservative. Therefore, the practitioner can use 
either (6) or the approximation (i). We recommend the use of approximation (i), as it is 
very simple to use and needs only the knowledge of standard normal percentile points. 

INSERT TABLE II HERE 

3. Estimation of reliability for repeated multivariate observations 

Let y x and x % denote the observation vectors (p- vectors) made on subject i at periods 1 
and 2 respectively. y l and Xi could also denote the observations on the same patients using 
two different treatments. Thus, y x and x, could be correlated. Assume that y ~ N\ (pi , Eu) 
and x ~ Ai(p 2 , E 2 ) and y l and x t may be correlated. Therefore, reliability is given by 

R 2 = P[y > x + r] = P[y - x > r], (10) 

Here x and y represents the p-variate random vector and r is a p-variate known quantity. 
The event [ y > x] represents that each component of y (i.e y t ) is greater than the corre- 
sponding component of x (i.e Xi). Without loss of generality we can assume that r = 0. 
Also, let d = y - x and d t — yi — x t , i — 1, . . . ,n, then d ~ N p (pd, E^), and R 2 = P[d > 0]. 
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Let p. d - d and t d = J ( d i ~ d){di ~ d ) T be the usual maximum likelihood estimators 

n i = i 

of p d and E d . Then, the MLE of R 2 is given by 


k? = p\d> o] w= ^ iEj= t„. (ii) 

It is known (see Anderson (1984)) that (d,S d ) are the complete sufficient statistics for 
(H d , E d ), where S d = ^E d . Now, let J (dl>0) be an indicator function such that I (dl>0) = 
1 , if ^ > 0 and = 0 otherwise. Then using an orthogonal transformation, Rao-Blackwell 
Theorem and Basu's Theorem, we show (see Appendix for the sketch of proof) that, the 
UMVUE of R 2 is given by 


£UMVUE = f f f f( D ) d D 

J V[di>0] ^ 

where 


( 12 ) 


D = 


J 


n 

n — 1 



(di ~ d ), 


f{D) = 


b 7 (1 -D T D)^f^ if 0 < D t D < 1 
0 otherwise, 


and 6 7 = T(^) [7r p / 2 • r(H^i)r . 

Obtaining UMVUE by (12) has to be done via numerical integration. A closed form 
estimate is not possible. However, since it has the desirable properties, namely, having the 
least variance among all the unbiased estimators, it is worth pursuing. We now conduct 
simulation studies to compare MLE and UMVUE of R 2 - 


Simulation Studies: 

We generated random samples of sizes n= 15,20,30, and 40, from bivariate normal 
distribution with several combinations of means, variances, and covariances. For each 
simulation, we calculated the MLE and UMVUE of R 2 , and their respective biases. At 
the end of 20,000 simulations, average bias, standard deviation and root mean square 
error for both the estimators are calculated. These are reported in Table III. As in the 
univariate case, once again the UMVUE seems to perform better than the MLE when the 
true probability is roughly between 0.25 and 0.75. However, there are a few cases when 
UMVUE seems to do better even when the true probability is about 0.11. It should also be 
noted that the difference in the mse of the two estimators is quite small once the sample size 
is large. Therefore, for large sample sizes MLE is preferred because of the simplicity in the 
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use. And for small to moderate sample sizes one has to choose either one based on whether 
one wants an unbiased estimate. We have not reported the actual values of the means, the 
standard deviations and correlations used in our simulations. Many combinations were 
chosen to obtain the true probabilities. Our main interest, however, is in studying the 
properties of MLE and UMVUE in relation to the true probability. Since, we are dealing 
with a bivariate normal distribution, it is clear that R 2 would be around 0.25 if the means 
are zeros and the correlation is zero (as the mass would be equally distributed in all four 
quadrants). Also, if the means are positive and the correlation is a high positive, much 
of the mass will be located in the first quadrant resulting in high true probability (i.e. 
.R 2 )- 1° practice, sample information available on means and correlations could be used in 
deciding whether to use MLE or UMVUE. Once again if one is interested in an estimator 
simple to calculate and that does quite well, MLE is the ideal choice. The computation 
of UMVUE for the bivariate case is quite simple, however, for higher dimensions it could 
take some time to calculate it, as it would involve integration of higher dimensions. With 
some programming knowledge one can calculate the UMVUE via subroutines from IMSL 
or NAG. 

INSERT TABLE III HERE 

4. Prediction interval for the future reliability when the observations on stress 
and strength are time dependent. 

Let (yi, x \ ),•••, (y n , x n ) denote the observations on stress and strength at time points 
1,2, and are time dependent (see Section 5.2 for examples). We assume that y 
and x are independent, and that the conditional distribution of (y n + i,£n+i) given y(n) = 
(yx, ■ ■ ■ , y n ) and x(n ) = (aq, • • • ,x n ) depends only on y n and x n . This assumption basically 
implies that there is a first order correlation (i.e. autoregression of order one) for both 
the processes { 2 } and {y}. We provide several examples (Section 5.2) to show that this 
assumption is satisfied easily. This assumption can be relaxed so that the conditional 
distribution depends on (y n -i,yn) and (z n _i, 2 n ), if the autoregression is of order two. 
Now, assuming that the observations have normal distribution, the future reliability at 
time (n + 1) given the past observations up to time n , is given by 

R 3 {yn,Xn) = P [yn +1 > ^n+i |y(n), rr(n)] 

[ioj 

= x ni9)\i 

where 9 is a vector of all the parameters, and 5 is some function of y n , x n and 9. 

Now, assume that there is enough information available, other than (y n ,£n), to esti- 
mate the parameter 9 and denote that estimator by 9 (see Kulkarni and Kushary (1991)). 
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This assumption is mostly required for the analytical derivations. We need this to obtain 
the conditional distribution of 5 — 8(y n , x n ,9) conditional on y n and x n , which is a crucial 
factor in constructing a lower bound on 6. However, in practice one can use y n and x n 
also in the estimation of 9 as they contribute highly to the prediction of y n + i and £n-fi- 
This should result in better estimates and herice yielding more precise results. 

Since reliability is conditional on y n and x n , we study the distribution of 8 — 8{y n , x n , 9) 
conditional on y n and x n . Now, using Taylor Series expansion of <5 as a function of 6, we 
have 


8(y n , x n , 9) = 8(y n ,Xn, 9) + {8'(y n ,x n , 9)) T [9 — 9) (14) 

where 8' denotes the derivative of 8 w.r.t. 9. 

Suppose there exists 9 such that 

(0-0)~iV(O,E n (6O). (15) 

Then, the conditional distribution of 8 given x n and y n is given by 

8(y n ,x n ,9)~N[8(y n ,x n ,9), al{9)} (16) 

where al{9) = [8'{y n , x n , 9)] T Z n (9)[8' (y n , x„, 9)}. Therefore, an (1 - a) 100% lower predic- 
tion bound for 8(y n ,x n ,9), given y n and x n , is 

L n = 8(y n ,x n , 9) - z a a n (9), (17) 

where a n {9) is a n {9) with 9 replaced by 9 and is a consistent estimator of a n {9). 

5. Applications 

5.1 Applications of procedures developed in Sections 2 and 3. 

We consider the data in Jordan et al. (1991), part of which is reproduced here in 
Table IV. They selected 10 critically ill patients that needed mechanical ventilation. Both 
the procedures (i.e. with NSI and without NSI) were administered on the same patients. 
Observations on heart rate, and ventilatory rate were made among other things. They 
recorded the measurements just prior to administering a procedure, right after adminis- 
tration and, 10 minutes after the procedure was administered. However, the measurements 
that are really of interest are prior to administration and right after administration, be- 
cause as more time is given patient’s conditions usually return to normal. After the first 
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administration, they allowed at least 90 minutes before administering the second proce- 
dure. The order of assignment was random. Another important variable measured was 

the amount of secretions removed. 

INSERT TABLE IV HERE 

First, we consider estimating the chances of 'with NSI being more effective in re- 
moving the secretions than 'without NSI' procedure. While the time allowed between the 
administration of the two procedures may have eliminated any carryover effects (from the 
previous procedure to the latter procedure), the patient effect is not removed. Therefore, 
the measurements made by the two procedures are paired observations. Let y and x denote 
the amount of secretions removed with and without NSI respectively. Then, for this data, 
R ™ = 0.99686 and rUMVUE _ q. 99988. Since, both MLE and UMVUE are very high, 
difference between them seems negligible. Therefore, one can choose either as their point 
estimate. Further, we construct a 95% most accurate lower confidence bound for Ri using 
all four methods. These are reported in Table V, along with the results for other two vari- 
ables of interest, namely, ventilatory rate and heart rate. However, for these variables, we 
define y =(post administration measurement with NSI — pre administration measurement 
with NSI), and similarly x is defined for the without NSI procedure. 

INSERT TABLE V HERE 

Next, since all the variables are measured simultaneously on each subject, they consti- 
tute a multivariable observation vector. Therefore, it is of interest to estimate the chances 
of patient being worse off with NSI compared to without NSI. It is known from the past 
studies, and clear in this study as well, that the amount of secretions removed using NSI is 
significantly more than without NSI. This is the main reason to administer NSI. Therefore, 
we will not include this variable in our multivariate study. The main purpose of the study 
is to see if there are increased chances of side effects (due to NSI) such as higher heart rate 
and ventilatory rates among other things. Therefore, we only consider two variables heart 
rate and ventilatory rate in our multivariate vector. Then, using the results of Section 3, 
estimates of reliability are given by R™= 0.41820, /jUMVUE = 0.40864. Therefore, there is 
a little over 40% chance that With NSI procedure results in elevated heart rate and ven- 
tilatory rate. In practice, just knowing the chance of increased heart rate and ventilatory 
rate may not be very informative as it does not say much about the magnitude of increase. 
If a physician determines that an increase of ‘a’ in the heart rate and an increase of ‘b’ 
in ventilatory rate might be dangerous, then one can calcuate the estimators according to 
that. Here we have simply calculated the chance of increased heart and ventilatory rate 
(i.e. a=b=0). 
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5.2 Examples for procedure in Section 4. 

First, we give several models that satisfy the conditional distribution assumption made 
in Section 4. 

Example 1: AR{ 1) model. 

Suppose y n and x n follow a stationary first order autoregressive model given by 


Vn — PlUn—l *b U n 


and 


x n = p 2 x n -i+v n (20) 

where u n and v n are independent normal r.v.’s with means 0 and variances a\ and o\. It 
is easy to see that the conditional distribution of (y n +i,x n+ i), given (y(n),x(n)) depends 
only on (y n iX n ) and of course the parameters. 

Example 2: Regression model with Autoregressive errors. 

Consider the usual regression model with errors following a first order autoregression, 
for both y and x, given by 


Pri — A)1 PllCin T Uni P\U n — \ "F Cln 


and 

X n = Po2 + ^12^2 n + v m V n = p2V n -l + C2n (21) 

where ei n and e 2n are independent normal random variables with means 0 and variances 
o\ and a 2 - This model also satisfies the condition that the conditional distribution of 
observations at time n + 1 depends only on the observations at time n. 

Example 3: First order autoregressive model with regressors. 

This model is similar to the model in Example 2, however one is not a simple 
reparametrization of the other (see Kulkarni (1987)). 


Vn = PlVn-l + @\\C\n + U n 


and 


X n — P2X n -l + Pl 2 C 2 n + Vm 


( 22 ) 


where u n and v n are independent normal random variables with means 0 and variances a\ 
and o\. It can be verified that assumption is satisfied for this model as well. 
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Next, we consider the model in (22) and illustrate the results of Section 4. Note, 
however, that similar results can be obtained for other models as well. For a good discussion 
of these models and related models see Basawa et al. (1985) and Kulkarni (1987). For 
simplicity, assume that o\ and o\ are known and let ci n = C2n • Then, 

A A A * 

<%n, x n j $) = Vn + Pl c n + 1 ” *2 x n ~~ I 

= 5{y n , x n , 0) + m T (9 - 9)/ 

where m T = (—y n ,—Cn+i>—Xn,—Cn+i), 9 T = (<5, fii, 62,02) = (0i\0j)> with 9j = (<5i ? A) 
and is the MLE of 9 obtained using ( 1 / 1 , Xi), • • • , (y n -i,Xn-i)- Notice that in (23) we did 
not have to use the Taylor series expansion as we could simply rewrite it in that form. If 
civ, + c 2n , then we would have to use the Taylor Series expansion. Also, assumption of o\ 
and o\ being known is for simplicity. When they are unknown, we incorporate these also in 
9 and obtain the asymptotic distribution of y/n- 1(0-0), as only (n- 1) observations are 
used in the estimation (see Kulkarni and Kushary (1991)). It can be shown using results 
in Kulkarni (1987) that (0* - 0)~N( 0, * = 1,2, where 

E* n (0i) = ^ 

cr i} 12 acr* 

n — 1 n — 1 n — 1 

with d = X) c ?i a i,i2 — yi-i c ij an( i a 2,i2 — ^37 XI ^-l^- Si = i(0t) is 

i= 1 7=1 ^ 7=1 

estimated consistently by replacing by Sf. Further, (5 — 5)~AT(<5, E n (0)), where 



E n (0) 


Eln(0) 0 

0 E 2n (0) ■ 


Using these results the conditional distribution of S(y n ,x n ,5), conditional on (y n ,Xn), Is 


S(y n , x n , $)~N(6(y n ,x n , S), <r*(0)), 

where a £(0) = ^-^rn T T, n (9)rn and m is as in (23). Also, estimated variance of S(y n ,x n ,9 n ) 
is given by a^(9). As noted earlier, in practice one can use y n and x n also in estimating 9 
and replace (n — 1) in all the expressions above by n. This should result in more precise 
results. 

At this point we should note that, even though a medical data set that would fit into 
our models is not presented here, the models described above are fairly common in medical 
research (see Diggle (1992)) and have broad applications. 
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Appendix: 

A sketch of the proof of the UMVUE for the p-variate normal case: 

Let xi ~ N p {fX\,Hx) and x 2 ~ N p (n 2 ,Y 2 ) and they may be are correlated. Consider, 
p = p(xi - x 2 > 0) = P(y > 0), where y = X\ - x 2 . Now onword, we will concentrate on 

the random variable y only. Let y ~ N p (y, X) and let y = y X = A Y (y» — y){yi — y) T 

i— 1 

denote the maximum likelihood estimators. Hence MLE of p is p = P(y > 0)|^_£ S _£. 
Let 6{yi) = 1, if y\ > 0 and 0 otherwise. Since (y, S ) are the complete sufficient statistics, 
where S — using Rao-Blackwell theorem E(6(yi)\y, S) is the UMVUE. Therefore, 

we need to find the distribution of y\ given (y, S). 

Joint distribution of j/i , • • • , y n is 

f(yi,--‘,yn) = Const • |X| ^ — tr|x 1 ^^{yi — n){yi — y-) 

n 

Consider the joint distribution of Z\, • ■ • , z n where Zi—Y^ Cjjyj, and Cj = ( Cn , • • • , Ci n ) 

j =1 

is the ith row of an orthogonal matrix whose with first and second rows are given by 
Cj = -A~\ T and Cj = c(l — £, — £, • ■ ■ respectively. Here, c is a normalizing constant 
given by c = (1 - A )" i. Using Theorem 3.3.1 of Anderson (1984), it can be shown that 

n 7i 

z\, z n are independent. Since, Y {Vi ~ p)(yi ~ P) T — Z) z i z T + { z i ~ V^y){zi — \/ny) , 

i— 1 i— 2 

77 n 

we have Y (Vi ~ V){Vi ~ vV — Z) Z * Z T ■ Thus, 

i= 1 i=2 


16 



s = y^ ZjZj + Z 2 Z% 

i= 3 

= 53 ZizJ + c 2 (yi - y)(yi - y) T 

i=3 

= S- 1 + c 2 (yi - y)(yi - y) T - 
Now the joint distribution of z\, ■ ■ • ,z n is 

rir n 

f(zi,Z 2 ,---z n )= Const • exp < --tr E -1 5~1 z i z I + ( z i ~ \/rcp)(ziy/n/i) T 

^ ^ L S=2 

Hence, zi, z 2 , S-i are independent with z\ ~ N(y/np,Yi) , z 2 ~ -N(0, E), and S-i ~ 
W( n _ 2,E), where W denotes the Wishart distribution. Now, the joint distribution of 

S- 1 and Z 2 is 

/(5 , _i,z 2 ) = Const •|E|"'^|S-i| (n_p_3)/2 exp{-±tr(E- 1 S_i)}|£| _ 5 
exp|-itr[E _l Z2^] 

From this, we find the joint distribution of S = SCi + 22 an< ^ z 2 given by 

f(S, z 2 ) = Const • |E|- |S - s 2 4'|<-'’- 3 > /2 exp|-itr[E-'S1 } 

where |5 — z 2 ^ 2 '| 1S positive definite. 

Further, the joint distribution of S and D = S~^Z 2 is 

r 1 (n— p-3)/ 2 r 1 ^ 

f(S,D) = Const-IEI -1 ^ \S\(I - DD T ) exp|--tr(E _1 5)| 

= Const ■ lEr^lSI 11 ^ ex p{-iir(E _1 5)J(l - D T D) {n - p ~ 3)/2 , 
implying that D is also independent of S. It can be seen that, 

f(D) = k(l- D T D) {n - p ~ 3)/2 , 0 < D t D < 1, 

where k = r( g ^)[7r p / 2 r( n ~^~ 1 )] -1 . 

Therefore, the distribution of D is free of the parameters. Since, the UMVUE of p is 
p = p(y 1 > Ojy, S) but the distribution of S~^(yi - y) is independent of parameters, by 
Basu’s theorem, conditional distribution of S~?(yi - y ) is independent of y and S. 

Hence, 

fi = j I f(D)dD. 

Vl>0 
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Table 1. Simulations results on mle (M) and UMVUE (U) with n=10 


True Probability 


Estimate 

Bias 

Std. dev. 

Root mse 

0.01222 

M 

0.01492 

0.00270 

0.02037 

0.02054 


U 

0.01225 

0.00003 

0.02190 

0.02190 

0.02275 

M 

0.02454 

0.00179 

0.02847 

0.02853 


U 

0.02241 

-0.00034 

0.03168 

0.03168 

0.04006 

M 

0.04088 

0.00082 

0.04071 

0.04072 


U 

0.04052 

0.00046 

0.04593 

0.04593 

0.06681 

M 

0.06394 

-0.00287 

0.05499 

0.05506 


U 

0.06635 

-0.00046 

0.06165 

0.06165 

0.10565 

M 

0.10006 

-0.00559 

0.07221 

0.07243 


U 

0.10636 

0.00071 

0.07898 

0.07898 

0.15866 

M 

0.14919 

-0.00946 

0.09077 

0.09126 


U 

0.15910 

0.00045 

0.09581 

0.09581 

0.22663 

M 

0.21463 

-0.01199 

0.10911 

0.10977 


U 

0.22655 

-0.00007 

0.11047 

0.11047 

0.30854 

M 

0.29653 

-0.01201 

0.12424 

0.12482 j 


U 

0.30759 

-0.00095 

0.12086 

0.12087 

0.40129 

M 

0.39639 

-0.00490 

0.13521 

0.13530 


U 

0.40294 

0.00164 

0.12775 

0.12776 

0.50000 

M 

0.50096 

0.00096 

0.13893 

0.13893 


U 

0.50093 

0.00093 

0.12989 

0.12989 

0.59871 

M 

0.60419 

0.00548 

0.13550 

0.13562 


U 

0.59760 

-0.00111 

0.12807 

0.12807 

0.69146 

M 

0.70348 

0.01202 

0.12382 

0.12440 


U 

0.69242 

0.00096 

0.12049 

0.12049 

0.77337 

Q 

0.78529 

0.01192 

0.10935 

0.11000 


H 

0.77337 

0.00000 

0.11064 

0.11064 

0.84134 

■MB 

0.85218 

EB 

0.09047 

0.09112 


H 

0.84238 

Ml 

0.09564 

0.09565 

0.89435 

M 

0.90059 

■ 

0.07172 

0.07199 


U 

0.89435 


0.07857 

0.07857 

0.93319 

M 

0.93613 

0.00293 

0.05462 

0.05470 


U 

0.93370 

0.00051 

0.06121 

0.06122 

0.95994 

M 

0.95970 

-0.00024 

0.04000 

0.04000 


U 

0.96014 

0.00020 

0.04512 

0.04512 

0.97725 

M 

0.97515 

-0.00210 

0.02903 

0.02911 


U 

0.97721 

-0.00004 

0.03233 

0.03233 

0.98778 

M 

0.98508 

-0.00269 

0.02026 

0.02044 


U 

0.98778 

0.00001 

0.02179 

0.02179 

0.99379 

M 

0.99133 

-0.00246 

0.01349 

0.01371 


U 

0.99391 

0.00012 

0.01362 

0.01362 







































































Table 2. Simulated coverages with true confidence 0.95. 


True Probability 

Using (6) 

Aprx.(i) 

Aprx.(ii) 

Aprx.(iii) 

n=10 





0.00620 

0.95000 

0.94745 

1.00000 

0.01855 

0.01222 

0.94810 

0.94425 

1.00000 

0.07605 

0.02275 

0.95190 

0.94965 

1.00000 

0.24055 

0.04006 

0.94960 

0.94745 

1.00000 

0.50180 

0.06681 

0.94965 

0.94850 

1.00000 

0.73925 

0.10565 

0.95075 

0.95000 

0.99985 

0.87680 

0.15866 

0.94980 

0.94955 

0.99850 

0.93240 

0.22633 

0.95215 

0.95215 

0.99355 

0.95030 

0.30854 

0.94880 

0.94880 

0.98155 

0.94845 

0.40129 

0.95065 

0.95070 

0.97485 

0.94560 

0.50000 

0.94780 

0.94830 

0.99150 

0.94245 

0.59871 

0.94860 

0.94975 

0.99915 

0.95055 

0.69146 

0.94930 

0.95105 

0.99990 

0.96835 

0.77337 

0.94970 

0.95295 

0.99995 

0.98215 

0.84134 

0.95440 

0.95790 

1.00000 

0.99540 

0.89435 

0.94920 

0.95375 

1.00000 

0.99820 

0.93319 

0.95005 

0.95485 

1.00000 

0.99970 

0.95994 

0.95435 

0.95980 

1.00000 

0.99995 

0.97725 

0.94885 

0.95685 

1.00000 

1.00000 

0.98778 

0.94880 

0.95625 

1.00000 

1.00000 

0.99379 

0.94630 

0.95285 

1.00000 

1.00000 







Table 3. Simulation Results for Bivariate Normal with n-15 


True Probability 

Estimators 

Bias 

St. Dev 

Root MSE 

0.0520 

M 

-0.0026 

0.0353 

0.0354 


U 

-0.0002 

0.0376 

0.0376 

0.1024 

M 

0.0002 

0.0446 

0.0446 


U 

0.0002 

0.0433 

0.0433 

0.1586 

M 

-0.0069 

0.0748 

0.0751 


U 

0.0002 

0.0776 

0.0776 

0.2015 

M 

0.0003 

0.0769 

0.0769 


U 

0.0002 

0.0743 

0.0743 

0.2500 

M 

-0.0008 

0.0893 

0.0893 


U 

-0.0008 

0.0861 

0.0861 

0.3500 

M 

0.0071 

0.0969 

0.0972 


U 

-0.0006 

0.0930 

0.0930 

0.4421 

M 

0.0116 

0.1003 

0.1010 


U 

-0.0013 

0.0974 

0.0974 

0.5075 

M 

0.0123 

0.1053 

0.1060 


U 

-0.0003 

0.1027 

0.1027 

0.5442 

M 

0.0126 

0.1061 

0.1068 


U 

0.0010 

0.1033 

0.1033 

0.6603 

M 

0.0075 

0.1004 

0.1007 


U 

-0.0011 

0.0988 

0.0988 

0.7079 

M 

0.0138 

0.0904 

0.0915 


U 

0.0001 

0.0933 

0.0933 

0.7525 

M 

0.0121 

0.0849 

0.0857 


U 

0.0003 

0.0887 

0.0887 

0.8105 

M 

0.0085 

0.0797 

0.0802 


U 

-0.0001 



0.8674 

M 

0.0039 




U 

-0.0003 


1 

0.9048 

M 


0.0512 

0.0513 


U 


0.0559 

0.0560 










Table 4. Data from Jordan and Garrett (1992) 


Patient 

Variable 

Time 

With NSI 

Without NSI 

i 

i 

Heart rate 

Pre 


122 



Post 


126 


Vent, rate 

Pre 

38 

42 



Post 

44 

50 


Secretions 

Post 

5 

0.5 

2 

Heart rate 

Pre 

81 

92 



Post 

111 

121 


Vent, rate 

Pre 

18 

14 



Post 




Secretions 

Post 



3 

Heart rate 

Pre 

120 

120 



Post 

125 

122 


Vent, rate 

Pre 

45 

50 



Post 

52 

60 


Secretions 

Post 

3.5 

2.0 

4 

Heart rate 

Pre 

84 

78 



Post 

89 

92 


Vent, rate 

Pre 

6 

8 



Post 

14 

12 


Secretions 

Post 

3.0 

1.5 

5 

Heart rate 

Pre 

132 

127 



Post 

144 

135 


Vent, rate 

Pre 

43 

42 



Post 

64 

56 


Secretions 

Post 

4.0 

2.0 

6 

Heart rate 

Pre 

121 

118 



Post 

140 

126 


Vent, rate 

Pre 

7 

8 

i 


Post 

14 

9 

| 

Secretions 

Post 

6 

4.0 

7 

Heart rate 

Pre 

135 

126 

| 


Post 

146 

130 


Vent, rate 

Pre 

43 

42 



Post 

64 

56 


Secretions 

Post 

3.5 

1.0 


8 

Heart rate 
Vent, rate 
Secretions 

9 

Heart rate 


Vent, rate 


Secretions 

10 

Heart rate 


Vent, rate 


Secretions 











































Table 5. Estimates and lower bounds for the data. 


Variable 

MLE 

UMVUE 

L-Using (6) 

Lai 



Secretions 
Heart Rate 
Vent. Rate 

0.99686 

0.62184 

0.57866 

0.99988 

0.61194 

0.57204 

1.43813 

-0.24557 

-0.34189 

1.40315 

-0.24582 

-0.34197 

0.37364 

-0.33987 

-0.46876 

0.96930 

-0.23241 

-0.33368 









